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Abstract 

In this paper is described a general 2"*^ order accurate (weak sense) pro- 
cedure for stablizing Monte-Carlo simulations of Ito stochastic differential 
equations. The splitting procedure includes explicit Runge-Kutta methods 
[1], semi-implicit methods [2] [3], and trapezoidal rule [1][4]. We prove the 
semi-implicit method of Ottinger [3] and note that it may be generalized for 
arbitrary splittings. 



1 Introduction 

^ We are interested in numerical procedures for simulating long time integrations of 

^ large dimensional Ito stochastic processes. The so-called Langevin dynamics ap- 

o proach of estimating physical quantities < / > uses a vector of stochastic processes 

^ a = 1, . . . ,n} which vary with simulated time t and converge to a station- 

ary state (e.g. Parisi and Wu [5], Klauder [6]). Estimates for < / > are long-time 
averages (large T) 



o 



<f>^^[dtf{x{t)), (1) 



00 

o 
o 
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^ where process x{t) satisfies a Langevin equation (Ito stochastic differential equa- 

0> tion) of multiplicative noise form 

dx" = r{x)dt + a"^(x) duj^{t), (2) 

driven by an n-dimensional vector of independent Brownian motions uj. Our no- 
tation in (2) uses the summation convention wherein repeated indices, in tliis case 
the /3 = 1, . . . ,n, are always summed over. Additionally, the following shorthand 
notation for partial derivatives will be found convenient: d^g = dg/dx'^ (for some 
function g{x)). In order that the long-time accuracy of (1) be of order of the time 
step size h, we desire that the single time step behavior of simulations of x{t) be 
of 0(/i2). 

Increments of the n independent Brownian motions u^{t) satisfy the relations 
{a, /3 = 1,... , n, where n is the size of vectors x and co, and S{t) is the Dirac 
delta function): 
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(j"(0) = (3) 

< duj^'it) > = 

< duj^{t2) > = (5"^ S{ti - 1,2) dtidh 

< du:^{t2) duj^{ts) dco^U) > = S^'^S^^ S{ti - ^2) S{t3 - U) dtidt2dhdU 

+ (5""'' 5{ti - h) S{t2 - U) dtidt2dtj.dti 

+ 5{ti - U) 5{t2 - ta) dtidt2dtsdU 

where the last relation follows from the third since the processes are Gaussian. 
These will suffice for all our needs to 0{h?'). 

The stochastic integral over a{x) depending on x must be handled carefully since 
the Brownian motion (Wiener process) ijj{t) is not of bounded variation. In this 
paper we choose the Ito definition wherein (2) is shorthand for 

= a;"(0) + f U'{x{s))ds + f a''^{x{s))dio^{s), (4) 

and the stochastic integral is belated, an Ito martingale [12]. In simple terms, the 
belated integral may be considered a Riemann sum in which the value taken for 
integrand a'^^{x{s)) in each t-interval {tj < s < tj+i) is for argument x{tj) at the 
beginning of the interval [6]. Other definitions, e.g. Stratonovich, are equivalent 
by transformations of the drift J bds under sufficient smoothness conditions on a 
(e.g. [12]). 



2 Numerical Approximations 

We begin by stating our algorithm and put off discussions of variants until later. 
For (2), writing b{x) split into two parts 

6"(x) = .4"(x) + 5"(a;), 

the following is a second order accurate (weak sense) simulation method. Vector 
xii is the process value at the end of time step h, and Xq is the process value at the 
beginning of the step. 

< +^(/l"(x;J+5"(xo + ao6 + (^o + 5o)/i) + ^"(:i:o) + 5"(xo)) 

+ \{a-'{xo + ^a,^o + ^{Ao + B,)) (5) 

+ ^"^(:^o-yfao6 + ^(A + 5o))Kf 
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In this formula (5), numbers .^o and .^i are vectors of independent, identically dis- 
tributed (iid), Gaussian (to 0{h?')) random variables with zero mean and variance 
h: 

<^o> = <er> = o, (6) 

<C?f> = 0, 

<^t^kUt> = {S'^^S^^ + S'''-'S'^^ + S''^S'^^)h\ i = 0,l. 

Again, the last relation follows from from the third since the .^o and are Gaussian. 
The are well modeled by = Vhzf", where zf are 2-n (q = 1, . . . , n and i = 0,1) 
normal random variates, each with zero mean and unit variance [1],[8],[9]. Array 
S"^^ contains models for the stochastic integrals (see Talay [11] for example, or 
Kloeden and Platen [6]): 

rt+h 

E'^ ^ uj'{s)dJ{s) = 0{h). 

These will be discussed in more detail in our proof of (5). Additionally, subscript 
on vector xq refers to the value of the process at the beginning of the time step 
t ^ t + h, and on the drift /diffusion coefficients indicate = A{xo), Bq = B{xo), 
etc. 

2.1 Taylor series and stochastic integrals 

In this section we set out some known results (e.g. Milstein book [9]) for the 
purposes of explaining notation and to make the arguments coherent. The idea of 
weak numerical approximations to (2) involves constructing a sequence of discrete 
xi.h 1 = 0,... vectors, one for each time step. Then, if xi.h = x{t), for an arbitrary 
smooth function / (C^ in each x"), we want the functionals 

< f{x{t + h))> = < f{x^i+iyh) > 

V ' ^ V ' 

computed using (3) computed using (6) 

to agree to some desired order in h. More precisely, if at time t, process x{t) has 
value xq and f{x{t)) has value /o, we require that 
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<f{x{t + hj) > 



+ 



2 dx'^dxi^ 




(7) 



+ 



1 



3! dx'^dxi^dx^ 



+ 



1 d'fo 



< Ax^Ax^Ax^Ax^ > 



4! dx^'dxi^dx^dx^ 



< f{xt+h) > +o{h'^). 



We compute the moments < Ax" • • • Ax^ > from expectation values of products 
of some stochastic integrals found when one step of (4) is expanded in a Taylor 
series. These moments must agree to the desired order in h {0{h?') in our case) 
with moments computed using the simulation (5). This idea, apparently due to 
Wagner and Platen [15], is well described in Milstein [9] or Kloeden and Platen [8]. 
That is, we write 



where AM"(/i) = 0{h^) + 0{h) + • • • , and A^"(/i) = 0{h) + 0(/J) + • • • , and 
x(s) appearing on the right hand side of (8) may be repeatedly substituted (in 
Picard fashion) by the right hand side of the integral formula (4). The result is 
a stochastic Taylor series (e.g. see [8], Chapter 5). If we can construct discrete 
models for the resultant stochastic integrals to 0(/i^), these models may be inserted 
into the stochastic Taylor series for Ax" to obtain a model Taylor series whose 
increments satisfy (7). 

2.2 Stochastic Taylor series 

We first derive the stochastic Taylor series to 0(/i^), then find some models for the 
required stochastic integrals. These results may be found elsewhere in more general 
formulations (e.g. [8]). It is important to note that x"(t) is a continuous Markov 
process and that our simulation is a discrete one. A method for computing any 
one step enables us to compute any other step. Hence, without loss of generality, 
the integrals from t ^ t + h my be abbreviated to ^ h. For example, consider 
the following step of the stochastic Taylor series: 





A^"(/i) 



AM"(/i) 
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/ a{x{s))duj{s) = / a{x{t) + Ax) d(jj{s) 
Jt Jt 

dco{s) + (9xO-o)cro Ij^ dco{u)> dco{s) + ... 

rt+h 

= (7o Aojh + {dx(Jo)(7o AujsdLo{s) + ... 

Now Acos = Jl^'duj{u) is a finite increment of Brownian motion co and has initial 
value Aloo = 0. Furthermore, the infinitesimal increments (in s) of AcOg are du:{s), 
whose properties are (3). Thus, Aa;^ may be treated as a Brownian motion in the 
finite interval t<t + s<t + h, i.e. < s < h, of interest. 

The two terms of (8) are 

AA'^ih) = b^h+{dBb^)a^^ dsco^s) (9) 

^0 

" V ' 

.r 

+ {dpb^) {dy,') a^' r ds r J {u) d^^iu) 



JO 



+ 



and 



AWih) = t d^^ih) +{9b(Jo'') (jt t ^"{s) duj^s) (10) 
Jq Jo 



fh fS 

H\ ::: :\ ::. : . 

=0 Ju=0 



Js=0 Ju=0 



+ 



Jo 



+ \ {dM'') (T^^ < [ ^"^^ duP + 0(^ 



In these expressions, those terms underbrace marked with * are 0{Y?) with van- 
ishing expectation values and may be ignored. That is, in (7) those 0{lr?) terms 
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whose expectation values are zero cannot contribute to any of the < Ax ■ ■ ■ Ax > 
moments to 0{h?'). In particular, all 0{h'^) terms of the martingale A^A"{h) may 
be ignored. The remaining labeled stochastic integrals are modeled as follows. 

2.3 Models for Stochastic Integrals 

In increasing order of powers of h (the step size), the needed stochastic integrals 
and their models are below. Following the enumeration, we prove the models only 
for and L'^'^^. Model S'^'^ is shown in [8] and [11] and like the remainder is a 
straightforward application of the correlations (3) and showing that products of the 
models have the same expectation values as their corresponding stochastic integrals 
to 0(/i^) accuracy. To verify the model for E"^^ , only products and times E"^^ , 
and S"^ times S'^'^ need be considered. No other terms of Ax in (8) can form any 
contribution to (7) to Oih^). 

0{h^): 

7" = Auj? = dco^is) ^ = Vhz^ 
Jo 

0{h'): 



Jo 



Here, variables the same zero mean, unit variance Gaussians 

that appear in the model for (respectively ^7)- Array z"^^ is defined 
only for e > 7 and consists oin- {n — l)/2 independent, zero mean, unit 
variance normal (to 0((5)^)) random variables. These are independent 
of the Zi, Zq, and each other. 
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0{hi): 



1 

2 

7 



Jo 

= [ sdio^(s) ^ IhCi 
Jo 2 

pn = co''{u)duj'{u))dco^{s) ^ 

Js=0 Ju=0 



or ^ UoeoCi 



In L'^'^^, the variables = ^/H zq contain independent zero mean, unit 
variance Gaussians z^: independent of tlie zf appearing in and the 
F'^ appearing in E"^^. 

['dsu;^{s)u;^{s) ^ 
Jo 2 

2 ' 

All three models for K'^'^ satisfy the calculus to 

Proof of model for J^: 

Since ~ |.^7 is 0{hi) with vanishing expectation, we observe that 
only products of .P with terms of 0{h^) can contribute to 0{h'^) in (7). 
Respectively then, 

• both P and the model have vanishing expectation, so the 
calculus is satisfied to 0{h^). 

• The expectation of product P times /" is 

<rr> = <{t dco"{u)){['' dsiT du'{v)))> 

Ju=0 Js=0 Jv=0 

rh rh rs 

= I du I ds I dvS°'^S{u — v) 

Ju=0 Js=0 Jv=0 

= -^"^ 

2 
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while the product of the models and | ^7 has expectation (in 
the space of z's): 

h h? h'^ 

<er^e7 >=y <^^7 >=y^- 

where (6) has been used. 

Proof of model for L'^'^'^: 

Since L"^'^'^ and the its proposed models are 0{h'^), we observe that 

• both L'^^'^ and the models have vanishing expectations. These 
follow from the facts that L'^'^"' is a martingale, and that the models 
are odd in in the Gaussian variables ^i, respectively. Thus, to 
0(/i2), the calculus is satisfied. 

• Inserting 

Lo^it) = f dco^iu), and 

Jv=Q 

into the expectation of product L'^'^^ times 7" we get 

Js=0 Jt=0 

/h rh rt rt 

ds dt dv du 
_=0 Jt=0 Jv=0 Ju=0 

{S'"'5'^5{s - u)5{t - v) 
+ - v)S{t - u) 

+ S'''^6'''6{s-t)6{u-v)} 

2 

In the four-fold integral: since v < t, the first term vanishes, 
and likewise u < t eliminates the second, so only the last term 
survives. The product of the model ^" and the first model |5°^'^^7 
has expectation (in the space of 2;'s): 

h h'^ h"^ 

SI 2 ^1 2 ^ 2 

Likewise, for the second proposed model, 

where (6) has been used to compute the expectations in both vari- 
ants. Thus, the models for L'^'^^ satisfy the calculus to 0{h?'). 
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2.4 Simplified Taylor series 

Substituting the models for stochastic integrals of section 2.3 into the Taylor series 
of section 2.2, we get a model Taylor series 

< = (11) 
+ Kh + {d0b^)at^^Ci drift A^"(/i) 

+ (^f^i+{df)(T^l(^t^"^ diffusion AM'^(/i) 

3 Proof of splitting formula 

Using the model Taylor series (11), we now verify the splitting formula (5). This 
is straightforward, writing for brevity (5) as 

Ax" = ^ + 5"(Xe„;er) + ^"(^^o) + 5"(Xo)} + AM"(/l) (12) 

where Xeuier = Xq-\- hbo + ctq^i is the Euler estimate, and 

= + {dpA^)Ax^ + l{dpd^A^)Ax^Ax^ + o{h) 

is all we need to 0(/i^). Expanding the second term on the right hand side one 
more time 

{dfiADAx^ = {dpA^){h{A^,+B^) + AM''{h)} + o{h) 
+ o{h) 

The 0{h) term containing may be ignored since it has vanishing expectation. 
The whole expression containing (9/3ylo)Ax^ is already 0{h), so this term will be 
0{h'^) overall. Subscripts on x, x+ and the arguments 
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4 = 4 + \ll<'(o+liA'o+B^,) 

Because 

a^^{x+) + a^^ix.) = 2ap + 0(/i), 
by way of the explicit construction of a:+ and x^ , we have 

{d^AD^x^ = (d^A-,) + ^0^) + ^^o'Ci] + o{h). 

Including this expansion in (12), 

Ax" = ^{A^ + (d^A^)(A^,+B^,)h+{df,A^)a^,^a + {dpd,A^)a^,'an^^^^^ 

B'^ix.^ler) + A'^iXo) + B"(xo)} + AM"(/l) 

leaving only the B{xeuier) term left. This is 

B^'iXeuler) = 5° (-^O + boh + (7o6) 

where 60 = Aq + Bq. We now need to expand this to 0{h), 

B^{x,^i,r) = B^ + {dpB^){A^o+B^)h+{dpB^)a^^Ci + 



+ 



\idpd,B^)ata^o'C,a + oih) 



Or, altogether, 



Ax" = ^{A-, + {d^A-,){A^o+B^o)h+{d^A-oyo'Ci 



2 

35 ^^e 



+{d8d,A^oWa^oriei + 

^0 + (^i?o )(4 + B^)h + (^i??)af ^7 

+ ^(a^a,5o")ao^Vo'''^^^i^ + A^ixo) + B^{xo)] 

AM"(/i) 



{A^o + i^o^)/^ + (5/5(^0 + B^A^o + B^o)y 
m^ + B^j)at'^a + 
{d^d,{A-o+B^))ata^o'CiCi + 
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Again using b = A + B, we get (11), which proves the formula (5). 

4 Explicit and semi- implicit variants, stability 

The significance of the splitting formula (5) lies principally in the increased stability 
of simulations when A 0. There are a plethora of references to implicit methods 
in ordinary differential equations (ODE's), most of which really only treat linear 
stability. For example, Bulirsch and Stoer [13] or Gear [14] discuss a basic analysis 
of the linear ODE, dx{t) = Axdt, where A is a negative matrix with eigenvalues 
having large scale differences. Namely, when |Amax|/|-^mm| is egregiously large, 
these scale differences in the eigenvalues (A) of A make simulations difficult. A 
time step h small enough to resolve short time scales is much too small to be 
practical for long time scale components. Increasing the time step leads to unstable 
simulations. Indeed, this basic linear analysis may be extended to the stochastic 
differential equation case. The obvious analog is a vector version of the Ornstein- 
Uhlenbeck process: dx = Axdt + duj{t). Several analysis of this situation exist (e.g. 
[1], [10]), and will be only briefiy discussed here. Instead, we will discuss some 
simple non-linear problems, where the drift is monomial. 

4.1 Trapezoidal rule and semi-implicit methods 

The solution of (5) can be a complicated affair when the drift is highly non-linear. 
At every time step, this equation must be solved for solution x^ which appears on 
both sides of the algorithm. However, in a common case 

6"(x) = + ^"(x), 

matrix A forms a linear part of b, and g{x) is non-linear. An obvious choice of 
splitting makes (5) easy to solve: 

A"(x) = A'^^x^ and B"(x) = g'^ix). 

This solution is affected by moving the ^Axh term which then appears on the right 
hand side of (5) to the left, computing (1 — |A)^^, and solving a linear system 
where everything on the right hand side is explicit. Variations similar to this have 
been known for a long time in ODE simulations (e.g. see Butcher [16]). For this 
semi-implicit splitting, the linear stability properties discussed in, say Gear [14], 
are preserved, but the resulting implicit equations are easy to solve. Other variants 
include the following. 

The choice 

A"{x) = b"{x) and S"(.'z;) = 0, 
gives an implicit trapezoidal rule, while the alternative 
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and 



A"{x) = 0, 



is an explicit trapezoidal rule, and is a 2 order Runge-Kutta method. Schurz 
[4] has shown that in many cases, in particular when the drift b = b{t) is not 
autonomous, some degree of imphcitness can be essential. His most glaring example 
is the case of Brownian bridge, x = io{t) — -^uj{ti), where x{0) = and x{ti) = 0, 
and to < t < ti. For that case, when A = ab{t) and B = {1 — a)b{t), only non-zero 
a (in fact, any non-zero a) gives correct results. Additionally, he showed that only 
implicit trapezoidal rule is asymptotically un-biased. 

Frequently, for example in polymer physics, it is the non-linear part of b that 
causes instabilities. Ottinger [3] has used a splitting wherein a linear part of b{x) = 
Ax + g{x), B{x) = Ax, is chosen for the explicit part, and the non-linear part, g{x), 
is taken for the implicit term A{x) = g{x). In the example in the next section (4.2), 
a monomial drift is illustrated wherein the increased stability is demonstrated. 

4.2 Simple stability analysis for implicit algorithms 

In simulations of ordinary differential equations, a method is said to be stable if 
when applied to the linear equation x = Ax, the discrete solution x{k ■ h) ^ 
as /c ^ oo when matrix A < 0. For Langevin equations, if a{x = 0) ^ 0, the 
solution doesn't degenerate, that is, x doesn't vanish even though b is contracting. 
The analog to x — > for the Langevin case is convergence to a strictly stationary 
process (e.g. see Doob [17]). Thus, we consider an approximate analysis of the 
following scalar additive noise problem 



where b is contracting. That is, b{x) is skew in the sense that b{x) < when a; >> 
and b(x) > when x << 0. As t ^ oo, the distribution function for x satisfies the 
forward Kolmogorov equation and becomes time independent: 



where N is a normalization. Let \x\'^^ = < |xp >t^oo be the asymptotic mean 
square. A necessary condition for mean square stability is then 



dx = b{x)dt + dw{t) 




Thus, the stationary distribution function in this limit is {x > 0) 




(13) 
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(14) 



whenever |a;oP >> la^p^^. In the infinitesimal h Umit, this condition for b to be 
contracting when |a:o| is large is 

2 Re {xob{xo)) + 1 < 0. 

For the discrete simulation, however, (14) is step size and process size dependent, 
just as in the ODE case. To illustrate the situation, we compare explicit and 
implicit trapezoidal rule algorithms. These are, from (5), respectively 

Xh=Xo + ^{b{Xeuler)+b{^o))+^, (15) 

where x^uier is the Euler estimate (see Section 3), and 



Xh = xo + -{b{xh) + b{xo)) + t (16) 

The necessary condition (14) applied to the two methods yields an approximate, 
and as we hopefully demonstrate, qualitatively correct, analysis obtained by ex- 
panding (15) and (16) to 0(/i^). This inequality (14) for the explicit form (15) 
is 



< l^fel^ > -|a;o|^ = \xo + hbo + Y^bol^ (1"^) 

+—{{x, + hbo + —b'Mb'^ 

< 0. 

We can get a semi-implicit approximation for the implicit trapezoidal rule as fol- 
lows. First, we move all the Xh dependent terms to the left hand side 

- ^K^h) =xo + ^b{xo) + ^, 
to be expanded in a Taylor series in Ax = Xh — xq. We get 

Ax - -b'^Ax - ^b'^{Axf = hbo + ^ + o{h^). 
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Now notice that {Ax^ = (^)^ + o{h), whence 



XH = xo + {i-lb',)-'Lbo + ^Ke+^]. 



(18) 



Squaring this expression and taking expectations produces a semi-implicit approx- 
imation for the inequality (14): 



< > -\xo\^ = 2{1 - h'o)-\hbo 



bo)xo 



(19) 



[i-^b',)-^hb,+^b'^^\ + \i-h',rh 



< 0. 

When a local linearization of the drift b{x) is permissible, in particular the Ornstein- 
Uhlenbeck process b{x) = b'^x, these conditions reduce to 

ii + K + ^(6'o)^rixor + ii + < ixor, 

or 



\l + hb', + ^%) 



I \2l2 



< 1. 



And, for the implicit case: 



(20) 



2°0 



\xo\^ + {l-h'^,)-^h< \Xo\\ 



or 



1 -I- !^h' 
i -t- 2^0 



r, Or, 



< 1. 



(21) 



We notice that although 6q < (6 is contracting), whenever {hb^l is large enough 
the first inequality (20) fails. However, for 6q < the second (21) is satisfied for 
large step-sizes. Experiments described in [1] show these conclusions in more detail. 

For non-linear problems where no Lipshitz bound on b{x) exists, and therefore a 
local linearization says nothing about larger xq behavior, the analysis isn't quite 
so simple but seems to give the same conclusion. Namely, the implicit rule (16) is 
significantly more stable than the explicit one (15). To put a finer point on this we 
looked at some monomial drift problems 
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dx = -xlxl'^'^dt + dwit), (22) 

which for integer m > quickly become stationary and (13) is easily evaluated 

(N fti 2(m+T)™^)' '^^^ ^^^^ m = is the Has'minskii process [18], which as t ^ oo 
becomes an exponential distribution {p{x,oc) ^ exp{—'2\x\)); and m = 1 is the 
Ornstein-Uhlenbeck or classic Langevin process which is asymptotically Gaussian. 
Computing the left-hand sides of inequalities (17) and (19), it was easy to study 
the stability regions. Our monomial drift examples have symmetric distributions, 
so only > is needed for illustration. Figures 1,2 plot the ratio 



< 


Xh\ 






\xo\ 


2 



for various step sizes {h = 1/10, 1/100) and values m = 2, 3, 4. To lowest order in 
h and small \xo\, this ratio is (/ ~ (l-'^^oP + ^O/l-'^oP- Note that < Ixh]"^ > — |a;oP = 
{q — 1) ■ |a;op, so g > 1 for large |a;o| indicates a diffusive or growing process, hence 
is unstable. Respectively, Figure 1 shows this function for the explicit algorithm 

(15) and Figure 2 that of (16). We note that clearly the implicit method keeps this 
function less than unity, and thus remains stable for all 

2 \ 
{m-l)h) 

This is the size limit at which the 0{h) diffusion term in (17), |1 + |6qP/i, vanishes. 
The explicit method (15) gives < |.Tftp > greater than l^op when |a:o| is large 
enough, and therefore becomes unstable. 

Finally, in the quadratic (m = 2) and cubic (m = 3) cases, the implicit formula 

(16) was solved exactly for Xh- Expanding the quadratic and cubic solutions in 
a Taylor series (in ^) permited an independent comparison with the approximate 
form (19). Namely, the expectation < > was computed to the desired order 
in < ^^■'^ >= 0(/i^) needed to achieve reasonable accuracy (e.g. plotting accuracy). 
Such comparisons for the m = 2, 3 cases are also shown in Figure 2 (labeled X2.1 
for m = 2 and h = 0.1, X2.01 for m = 2 and h = 0.01, etc.). These X2.1, X2.01 
(m = 2), and X3.1, X3.01 (m = 3) exact solutions of the discretized equation (16) 
are to be compared with the 2.1, 2.01 and 3.1, 3.01 curves, respectively, which 
were computed using the approximate formula (19). We note that whenever |a;o| is 
not too large, the approximation (19) gives quite good results, and is qualitatively 
correct for larger process values. A fortiori, the exact solutions shown by curves 
X2.1, X2.01, X3.1, and X3.01 actually show better stability than the approximation 
(19) calculated using (18) shown by curves 2.1, 2.01, 3.1, and 3.01 in Figure 2. For 
the monomial drifts, thus polynomial drifts, estimate (19) thus gratefully appears 
pessimistic. 
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100 



Stability ratio vs. |xq| 

[dx = -x|xr"^dl+dw, m = 2,3,4] 




J I I I I I I I I I I I I L 

0.1 1 10 100 

process value Ixgl 



Figure 1: Stability of explicit trapezoidal rule (15) for problem (22). Curves labeled 
m.d mean m of (22) with stepsize h = O.d = 0.1, or 0.01. 
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100 



Stability ratio vs. IxqI 

[dx = -x|xr"Ml+dw, m = 2,3,4] 




0.1 1 10 100 

process value [xqI 



Figure 2: Stability of implicit trapezoidal rule (15) for problem (22). Curves labeled 
m.d mean m of (22) with stepsize h = O.d = 0.1 or 0.01. Labels Xm.d refer to 
exact solutions of the implicit difference equation (16) and are to be compared to 
the approximate analysis (19). 
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5 Conclusions 



We have shown a general splitting for a 2nd order weak accurate simulation method 
for Ito stochastic differential equations. The splitting permits choices for implicit 
dependences in the discretized time stepping which can improve stability. These 
choices can be made according to ease of solution (linear semi-implicit methods), 
or to improve stability when non-linear drift terms cause difficulties. Such meth- 
ods have been shown useful in polymer physics [3]. Additionally, an approximate 
analysis of stability computed using the semi-implicit approximation (18) seems to 
yield quantitatively reliable predictions for additive noise problems when process 
sizes aren't too large, and seems qualitatively reliable in any case. This analysis 
shows that improvement in stability can be expected at least for polynomial non- 
linearities by using implicit trapezoidal rule. Finally, it has not escaped our notice 
that (18) is easily generalizable to the multiplicative noise case {cr{x) 7^ 1), to yield 
a 2nd order weakly accurate linearly stabilized algorithm with no implicit equa- 
tions to solve. The drawback of such a procedure being principally more functional 
evaluations (i.e. b',b", if available), but will still be doubtlessly easier than solving 
implicit equations in the general case. 
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